
#This project is statistical analysis for a 3X1 factorial experiment-1 and 3x2 factorial experiment-2. For experiment-1 contain four scenarios for each of three conditions.
# The conditions involve AI, Human and Human aided by AI. each of these are dummy coded # to run the regression model.
# TRUST, CREDEBILITY, AGREEABLENESS, AND FAIRNESS are the dependent variables for experiment-1. 
# TRUST, CREDEBILITY, QUALITY, AND BIAS are the dependent variables for experiment-2 as perceived by #participants. conditions is independent variable.

# Step-0: pre-process data to remove participants who does not qualify
# Step-1: regression analysis for experiment-1 and 2
# Step-2: MLM models experiment-1.
# Step-2: MLM models experiment-2.

########################################################################
#                                                                      #
#                     0 Data cleaning process                          #
#                                                                      #
########################################################################
library(MASS)
library(dplyr)
library(scales)
library(performance)
library(reshape2)
library(nlme)
library(car)
library(MuMIn)
library(emmeans)


#reading the data file (appropriate file name goes there)
ex.data=read.csv("xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx")

#identifying the good completes who passed the filter questions and attention check
dat_com = ex.data %>% filter((SOCIAL_1!=1 | SOCIAL_2!=1 | SOCIAL_3 !=1 | SOCIAL_4!=1 ) & CHECK==2)
nrow(dat_com)


# computing the time quantiles to identify speeders
dat_com$Duration..in.seconds. <- as.numeric(dat_com$Duration..in.seconds.)
class(dat_com$Duration..in.seconds.)
dat_com_str<- dat_com
median(as.numeric(ex.data$Duration..in.seconds.))
quantile(ex.data$Duration..in.seconds.)


# excluding-speeders who took below 48% of median time
dat_com_str_sp<-dat_com_str %>%
  filter(Duration..in.seconds. > 268.60)
nrow(dat_com_str_sp)


# excluding duplicates
final_data <- dat_com_str_sp %>% distinct(pid, .keep_all = T)
nrow(final_data)

# excluding incomplete responses
final_data <- final_data %>% filter(!is.na(final_data$BEN_CRED))
nrow(final_data)



##################################################################################################
#                                                                                                #
#                   1 Regression analysis for experiment-1                                       #
#                                                                                                #
##################################################################################################
########################################################################
#                                                                      #
#               1.1 Creating variables of interest                     #
#                                                                      #
########################################################################
#### creating dummy variables (IVs) for AI, Humnan, and Mix conditions 
final_data$AI.dummy <- ifelse(final_data$condition %in% c(2,3),0,1)
final_data$Human.dummy <- ifelse(final_data$condition %in% c(1,3),0,1)
final_data$AI_Human.dummy <- ifelse(final_data$condition %in% c(1,2),0,1)
# Here AI condition is 1, Human condition is 2, and Mix condition is 3

#### Creating aggregated dependent variables (trust, fair, bias, and agree) across scenarios
final_data$mean.trust <- (final_data$TRUST+final_data$CIVIL_TRUST+ final_data$FACTS_TRUST+ final_data$DEC_TRUST)/4
final_data$mean.fair <- (final_data$FAIR+final_data$CIVIL_FAIR+ final_data$FACTS_FAIR+ final_data$DEC_FAIR)/4
final_data$mean.biased <- (final_data$BIASED+final_data$CIVIL_BIASED+ final_data$FACTS_BIASED+ final_data$DEC_BIASED)/4
final_data$mean.agree <- (final_data$AGREE+final_data$CIVIL_AGREE+ final_data$FACTS_AGREE+ final_data$DEC_AGREE)/4


### Creating aggregated dependent variables for each scenarios

#scenario1
final_data$agg_s1 <- (final_data$TRUST + final_data$FAIR + final_data$AGREE + final_data$BIASED)/4 

#scenario2
final_data$agg_s2 <- (final_data$CIVIL_TRUST + final_data$CIVIL_FAIR + final_data$CIVIL_AGREE + final_data$CIVIL_BIASED)/4 

#scenario3
final_data$agg_s3 <- (final_data$FACTS_TRUST + final_data$FACTS_FAIR + final_data$FACTS_AGREE + final_data$FACTS_BIASED)/4  

#scenario4
final_data$agg_s4 <- (final_data$DEC_TRUST + final_data$DEC_FAIR + final_data$DEC_AGREE + final_data$DEC_BIASED)/4 

### In addition to these variables, AI_awareness variable was hand coded as it was a multiple answers to be aggregated


########################################################################
#                                                                      #
#     1.2 Randomization check across different control variables       #
#                                                                      #
########################################################################

final_data %>% group_by(condition) %>% summarise(mean(RACE), sd(RACE))
final_data %>% group_by(condition) %>% summarise(mean(GENDER), sd(GENDER))
final_data %>% group_by(condition) %>% summarise(mean(EDU, na.rm=T), sd(EDU, na.rm=T))
final_data %>% group_by(condition) %>% summarise(mean(AGE_CAT, na.rm=T), sd(AGE_CAT, na.rm=T))
final_data %>% group_by(condition) %>% summarise(mean(PARTY7, na.rm=T), sd(PARTY7, na.rm=T))
final_data %>% group_by(condition) %>% summarise(mean(IDEO_1, na.rm=T), sd(IDEO_1, na.rm=T))
final_data %>% group_by(condition) %>% summarise(mean(FOLLOW, na.rm=T), sd(FOLLOW, na.rm=T))

#here AI is AI_awareness variable which is handcoded in the processed datasrt
#final_data %>% group_by(condition) %>% summarise(mean(AI, na.rm=T), sd(AI, na.rm=T))
#final_data %>% group_by(condition) %>% summarise(mean(as.numeric(AI_2), na.rm=T), sd(AI_2, na.rm=T))

summary(aov(scale(condition)~scale(RACE), data = subset(final_data, !is.na(RACE)) ), na.rm=T)
summary(aov(scale(condition)~scale(GENDER) , data = subset(final_data,!is.na(GENDER))), na.rm=T)
summary(aov(scale(condition)~scale(EDU), data = subset(final_data, !is.na(EDU))), na.rm=T)
summary(aov(scale(condition)~scale(AGE_CAT), data = subset(final_data, !is.na(AGE_CAT))), na.rm=T)
summary(aov(scale(condition)~scale(PARTY7), data = subset(final_data, !is.na(PARTY7))), na.rm=T)
summary(aov(scale(condition)~scale(IDEO_1), data = subset(final_data, !is.na(IDEO_1))), na.rm=T)
summary(aov(scale(condition)~scale(FOLLOW), data = subset(final_data, !is.na(FOLLOW))), na.rm=T)
summary(aov(scale(condition)~scale(INCOME), data = subset(final_data, !is.na(FOLLOW))), na.rm=T)
#summary(aov(scale(condition)~scale(AI), data = subset(final_data, !is.na(AI))), na.rm=T)

### randomization check revealed that three conditions are significantly differ for education variable. Hence we control for it in our models.


########################################################################
#                                                                      #
# 1.3 Regression analysis for study1 (aggregated dependent variables)  #
#                                                                      #
########################################################################

#human as reference conditions
summary(lm((mean.trust) ~ (AI.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((mean.fair) ~ (AI.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((mean.biased) ~ (AI.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((mean.agree) ~ (AI.dummy) + (AI_Human.dummy)+(EDU), data = final_data), na.rm=T)

# AI as reference conditions
summary(lm((mean.trust) ~ (Human.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((mean.fair) ~ (Human.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((mean.biased) ~ (Human.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((mean.agree) ~ (Human.dummy) + (AI_Human.dummy)+(EDU), data = final_data), na.rm=T)

#AI awareness as moderator
summary(lm((mean.trust) ~ (AI.dummy) + (AI_Human.dummy) +(EDU)+ AI + AI*AI.dummy + AI*AI_Human.dummy, data = final_data), na.rm=T)
summary(lm((mean.fair) ~ (AI.dummy) + (AI_Human.dummy) +(EDU)+ AI + AI*AI.dummy + AI*AI_Human.dummy, data = final_data), na.rm=T)
summary(lm((mean.biased) ~ (AI.dummy) + (AI_Human.dummy) +(EDU)+ AI + AI*AI.dummy + AI*AI_Human.dummy, data = final_data), na.rm=T)
summary(lm((mean.agree) ~ (AI.dummy) + (AI_Human.dummy)+(EDU)+ AI + AI*AI.dummy + AI*AI_Human.dummy, data = final_data), na.rm=T)

summary(lm((mean.trust) ~ (Human.dummy) + (AI_Human.dummy) +(EDU)+ AI + AI*Human.dummy + AI*AI_Human.dummy, data = final_data), na.rm=T)
summary(lm((mean.fair) ~ (Human.dummy) + (AI_Human.dummy) +(EDU)+ AI + AI*Human.dummy + AI*AI_Human.dummy, data = final_data), na.rm=T)
summary(lm((mean.biased) ~ (Human.dummy) + (AI_Human.dummy) +(EDU)+ AI + AI*Human.dummy + AI*AI_Human.dummy, data = final_data), na.rm=T)
summary(lm((mean.agree) ~ (Human.dummy) + (AI_Human.dummy)+(EDU)+ AI + AI*Human.dummy + AI*AI_Human.dummy, data = final_data), na.rm=T)


########################################################################
#                                                                      #
# 1.4 Regression analysis for study1 (Scenariowise for each dv)        #
#                                                                      #
########################################################################

######################################### Scenario-1 ######################################################################
summary(lm((TRUST) ~ (AI.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((FAIR) ~ (AI.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((BIASED) ~ (AI.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((AGREE) ~ (AI.dummy) + (AI_Human.dummy)+(EDU), data = final_data), na.rm=T)

summary(lm((TRUST) ~ (Human.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((FAIR) ~ (Human.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((BIASED) ~ (Human.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((AGREE) ~ (Human.dummy) + (AI_Human.dummy)+(EDU), data = final_data), na.rm=T)

#AI awareness as moderator
summary(lm((TRUST) ~ (AI.dummy) + (AI_Human.dummy) +(EDU)+ AI + AI*AI.dummy + AI*AI_Human.dummy, data = final_data), na.rm=T)
summary(lm((FAIR) ~ (AI.dummy) + (AI_Human.dummy) +(EDU)+ AI + AI*AI.dummy + AI*AI_Human.dummy, data = final_data), na.rm=T)
summary(lm((BIASED) ~ (AI.dummy) + (AI_Human.dummy) +(EDU)+ AI + AI*AI.dummy + AI*AI_Human.dummy, data = final_data), na.rm=T)
summary(lm((AGREE) ~ (AI.dummy) + (AI_Human.dummy)+(EDU)+ AI + AI*AI.dummy + AI*AI_Human.dummy, data = final_data), na.rm=T)

summary(lm((TRUST) ~ (Human.dummy) + (AI_Human.dummy) +(EDU)+ AI + AI*Human.dummy + AI*AI_Human.dummy, data = final_data), na.rm=T)
summary(lm((FAIR) ~ (Human.dummy) + (AI_Human.dummy) +(EDU)+ AI + AI*Human.dummy + AI*AI_Human.dummy, data = final_data), na.rm=T)
summary(lm((BIASED) ~ (Human.dummy) + (AI_Human.dummy) +(EDU)+ AI + AI*Human.dummy + AI*AI_Human.dummy, data = final_data), na.rm=T)
summary(lm((AGREE) ~ (Human.dummy) + (AI_Human.dummy)+(EDU)+ AI + AI*Human.dummy + AI*AI_Human.dummy, data = final_data), na.rm=T)

############################################ Scenario-2 ####################################################################
summary(lm((CIVIL_TRUST) ~ (AI.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((CIVIL_FAIR) ~ (AI.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((CIVIL_BIASED) ~ (AI.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((CIVIL_AGREE) ~ (AI.dummy) + (AI_Human.dummy)+(EDU), data = final_data), na.rm=T)

summary(lm((CIVIL_TRUST) ~ (Human.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((CIVIL_FAIR) ~ (Human.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((CIVIL_BIASED) ~ (Human.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((CIVIL_AGREE) ~ (Human.dummy) + (AI_Human.dummy)+(EDU), data = final_data), na.rm=T)

############################################## Scenario-3 #############################################################

summary(lm((FACTS_TRUST) ~ (AI.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((FACTS_FAIR) ~ (AI.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((FACTS_BIASED) ~ (AI.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((FACTS_AGREE) ~ (AI.dummy) + (AI_Human.dummy)+(EDU), data = final_data), na.rm=T)

summary(lm((FACTS_TRUST) ~ (Human.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((FACTS_FAIR) ~ (Human.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((FACTS_BIASED) ~ (Human.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((FACTS_AGREE) ~ (Human.dummy) + (AI_Human.dummy)+(EDU), data = final_data), na.rm=T)

################################################## Scenario-4 #########################################################

summary(lm((DEC_TRUST) ~ (AI.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((DEC_FAIR) ~ (AI.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((DEC_BIASED) ~ (AI.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((DEC_AGREE) ~ (AI.dummy) + (AI_Human.dummy)+(EDU), data = final_data), na.rm=T)

summary(lm((DEC_TRUST) ~ (Human.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((DEC_FAIR) ~ (Human.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((DEC_BIASED) ~ (Human.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((DEC_AGREE) ~ (Human.dummy) + (AI_Human.dummy)+(EDU), data = final_data), na.rm=T)


########################################################################
#                                                                      #
#    1.5 Regression analysis for study1 (aggregated DV)                #
#                                                                      #
########################################################################


## scenario1
summary(lm((agg_s1) ~ (AI.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((agg_s1) ~ (Human.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)

## scenario2
summary(lm((agg_s2) ~ (AI.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((agg_s2) ~ (Human.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)

## scenario3
summary(lm((agg_s3) ~ (AI.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((agg_s3) ~ (Human.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)

## scenario4
summary(lm((agg_s4) ~ (AI.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)
summary(lm((agg_s4) ~ (Human.dummy) + (AI_Human.dummy) +(EDU), data = final_data), na.rm=T)


########################################################################
#                                                                      #
#         1.6 Internal consistency check and factor analysis            #
#                             Appendix-5 and 7                         #
########################################################################
# calculating cronbach's alpha for internal consistency for US sample
# use the subset of the questions whose internal consistency is to be calculated

library(psych)

#1. taking subset of DVs for each scenario for study1
alpha_s1s1 <- final_data[, 49:52]
alpha(alpha_s1s1)

alpha_s1s2 <- final_data[, 53:56]
alpha(alpha_s1s2)

alpha_s1s3 <- final_data[, 57:60]
alpha(alpha_s1s3)

alpha_s1s4 <- final_data[, 61:64]
alpha(alpha_s1s4)

## taking subset of DVs for each issue for study2

alpha_s2s1 <- final_data[, 65:68]
alpha(alpha_s2s1)

alpha_s2s2 <- final_data[, 69:72]
alpha(alpha_s2s2)

alpha_s2s3 <- final_data[, 73:76]
alpha(alpha_s2s3)

# 2.Choose the "right" number of factors to retain
##parallel analysis
library(nFactors)
library(psych)
library(GPArotation)

correl = cor(alpha_s1s1, use = "pairwise.complete.obs")
correl
symnum(correl)

ev <- eigen(cor(alpha_s1s1)) # get eigenvalues 
ap <- parallel(subject=nrow(x = alpha_s1s1),var=ncol(x = alpha_s1s1), rep=100, cent=.05) 
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) 
plotnScree(nS)
ap$eigen
nofactors1 <- fa.parallel(alpha_s1s1, fm='ml', fa= 'fa')
sum(nofactors1$fa.values >.7)


s1.fa1 <- fa(r = cor(alpha_s1s1), nfactors = 2, rotate = "oblimin", fm = "ml")
plot(s1.fa1,labels=names(x = alpha_s1s1),cex=.7, ylim=c(-.1,1)) 
s1.fa1

############### scenario2############
correl = cor(alpha_s1s2, use = "pairwise.complete.obs")
symnum(correl)

# 2.Choose the "right" number of factors to retain
##parallel analysis

ev <- eigen(cor(alpha_s1s2)) # get eigenvalues 
ap <- parallel(subject=nrow(x = alpha_s1s2),var=ncol(x = alpha_s1s2), rep=100, cent=.05) 
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) 
plotnScree(nS)
ap$eigen

nofactors2 <- fa.parallel(alpha_s1s2, fm='ml', fa= 'fa')
sum(nofactors2$fa.values >.7)


s1.fa2 <- fa(r = cor(alpha_s1s2), nfactors = 2, rotate = "oblimin", fm = "ml")
plot(s1.fa2,labels=names(x = alpha_s1s2),cex=.7, ylim=c(-.1,1)) 
s1.fa2

############### scenario3 ############
correl = cor(alpha_s1s3, use = "pairwise.complete.obs")
symnum(correl)

# 2.Choose the "right" number of factors to retain
##parallel analysis
library(nFactors) 
ev <- eigen(cor(alpha_s1s3)) # get eigenvalues 
ap <- parallel(subject=nrow(x = alpha_s1s3),var=ncol(x = alpha_s1s3), rep=100, cent=.05) 
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) 
plotnScree(nS)
ap$eigen

nofactors3 <- fa.parallel(alpha_s1s3, fm='ml', fa= 'fa')
sum(nofactors3$fa.values >.7)


s1.fa3 <- fa(r = cor(alpha_s1s3), nfactors = 2, rotate = "oblimin", fm = "ml")
plot(s1.fa3,labels=names(x = alpha_s1s3),cex=.7, ylim=c(-.1,1)) 
s1.fa3


########################## Scenario4 #######################

correl = cor(alpha_s1s4, use = "pairwise.complete.obs")
symnum(correl)

# 2.Choose the "right" number of factors to retain
##parallel analysis

ev <- eigen(cor(alpha_s1s4))# get eigenvalues 
ev$values
ap <- parallel(subject=nrow(x = alpha_s1s4),var=ncol(x = alpha_s1s4), rep=100, cent=.05) 
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) 
plotnScree(nS)
ap$eigen

nofactors4 <- fa.parallel(alpha_s1s4, fm='ml', fa= 'fa')
sum(nofactors4$fa.values >.7)


s1.fa4 <- fa(r = cor(alpha_s1s4), nfactors = 2, rotate = "oblimin", fm = "ml")
plot(s1.fa4,labels=names(x = alpha_s1s4),cex=.7, ylim=c(-.1,1)) 
s1.fa4

#########################Study2- Issue1 ###################

correl = cor(alpha_s2s1, use = "pairwise.complete.obs")
symnum(correl)

# 2.Choose the "right" number of factors to retain
##parallel analysis

ev <- eigen(cor(alpha_s2s1)) # get eigenvalues 
ap <- parallel(subject=nrow(x = alpha_s2s1),var=ncol(x = alpha_s2s1), rep=100, cent=.05) 
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) 
plotnScree(nS)
ap$eigen

nofactors1 <- fa.parallel(alpha_s2s1, fm='ml', fa= 'fa')
sum(nofactors1$fa.values >.7)


s2.fa1 <- fa(r = cor(alpha_s2s1), nfactors = 2, rotate = "oblimin", fm = "ml")
plot(s2.fa1,labels=names(x = alpha_s2s1),cex=.7, ylim=c(-.1,1)) 
s2.fa1

#########################Study2- Issue2 ###################

correl = cor(alpha_s2s2, use = "pairwise.complete.obs")
symnum(correl)

# 2.Choose the "right" number of factors to retain
##parallel analysis

ev <- eigen(cor(alpha_s2s2)) # get eigenvalues 
ap <- parallel(subject=nrow(x = alpha_s2s2),var=ncol(x = alpha_s2s2), rep=100, cent=.05) 
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) 
plotnScree(nS)
ap$eigen

nofactors2 <- fa.parallel(alpha_s2s2, fm='ml', fa= 'fa')
sum(nofactors2$fa.values >.7)


s2.fa2 <- fa(r = cor(alpha_s2s2), nfactors = 2, rotate = "oblimin", fm = "ml")
plot(s2.fa2,labels=names(x = alpha_s2s2),cex=.7, ylim=c(-.1,1)) 
s2.fa2

#########################Study2- Issue3 ###################

correl = cor(alpha_s2s3, use = "pairwise.complete.obs")
symnum(correl)

# 2.Choose the "right" number of factors to retain
##parallel analysis

ev <- eigen(cor(alpha_s2s3)) # get eigenvalues 
ap <- parallel(subject=nrow(x = alpha_s2s3),var=ncol(x = alpha_s2s3), rep=100, cent=.05) 
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) 
plotnScree(nS)
ap$eigen

nofactors3 <- fa.parallel(alpha_s2s3, fm='ml', fa= 'fa')
sum(nofactors3$fa.values >.7)


s2.fa3 <- fa(r = cor(alpha_s2s3), nfactors = 2, rotate = "oblimin", fm = "ml")
plot(s2.fa3,labels=names(x = alpha_s2s3),cex=.7, ylim=c(-.1,1)) 
s2.fa3


#############################################################################################
#                                                                                           #
#                   2  Multilevel models for experiment-1                                   #
#                                                                                           #
#                                                                                           #
#############################################################################################
########################################################################
#                                                                      #
#               2.1 Creating variables of interest                     #
#                                                                      #
########################################################################



###1. Load the data and assign the unique ID to participants. Generate the aggregated DV for all scenarios.

#creating datset backup
mlm_data_us <- final_data

#assinging ID to each participant
mlm_data_us$ID <- 1:nrow(mlm_data_us)
mlm_data_us$agg_s1 <- (mlm_data_us$TRUST + mlm_data_us$FAIR + mlm_data_us$AGREE + mlm_data_us$BIASED)/4
mlm_data_us$agg_s2 <- (mlm_data_us$CIVIL_TRUST + mlm_data_us$CIVIL_FAIR + mlm_data_us$CIVIL_AGREE + mlm_data_us$CIVIL_BIASED)/4
mlm_data_us$agg_s3 <- (mlm_data_us$FACTS_TRUST + mlm_data_us$FACTS_FAIR + mlm_data_us$FACTS_AGREE + mlm_data_us$FACTS_BIASED)/4
mlm_data_us$agg_s4 <- (mlm_data_us$DEC_TRUST + mlm_data_us$DEC_FAIR + mlm_data_us$DEC_AGREE + mlm_data_us$DEC_BIASED)/4


### 2	Convert the wide data format to long data as MLM uses wide data format. Since, we have four scenarios we
###       consider each scenario as level. hence the data has four levels for study1.

mlm_data_us_L <- melt(mlm_data_us, id.vars = c("ID", "EDU"), measure.vars = c("agg_s1", "agg_s2", "agg_s3", "agg_s4"))
mlm_data_us_L2 <- melt(mlm_data_us, id.vars = c("ID", "condition", "EDU") , measure.vars = c("agg_s1", "agg_s2", "agg_s3", "agg_s4"))

mlm_data_us_L2$scenario <- recode(mlm_data_us_L2$variable, '"agg_s1" = 1; "agg_s2" = 2;
"agg_s3" = 3; "agg_s4" = 4;')

#Create dummy variables for conditions

mlm_data_us_L2$AI <- ifelse(mlm_data_us_L2$condition == 1, 1, 0)
mlm_data_us_L2$AI_Human <- ifelse(mlm_data_us_L2$condition == 3, 1, 0)
mlm_data_us_L2$Human <- ifelse(mlm_data_us_L2$condition == 2, 1, 0)
mlm_data_us_L2$moderate <- ifelse(mlm_data_us_L2$scenario == 1, 1, 0)
mlm_data_us_L2$civil_reminder <- ifelse(mlm_data_us_L2$scenario == 2, 1, 0)
mlm_data_us_L2$present_facts <- ifelse(mlm_data_us_L2$scenario == 3, 1, 0)
mlm_data_us_L2$news <- ifelse(mlm_data_us_L2$scenario == 4, 1, 0)

########################################################################
#                                                                      #
#           2.2 Multilevel models and main analysis                    #
#                                                                      #
########################################################################


###Multilevel Models for testing H1 and H2

##Human and news recommendation as the reference categories (H1, Appendix 8. Table 1) 

model_e1_us <- lme(value ~ AI + AI_Human + moderate + civil_reminder + present_facts +
                     AI*moderate + AI*civil_reminder + AI*present_facts + 
                     AI_Human*moderate + AI_Human*civil_reminder + AI_Human*present_facts + scale(EDU),
                   random =~ 1|ID, data=mlm_data_us_L2)
summary(model_e1_us)
anova(model_e1_us)
r.squaredGLMM(model_e1_us)

##Estimated marginal means (Appendix 8. Table 2)

#Recoding scenario instead of using dummy variables to make visualization easier (coding news as reference cat).
mlm_data_us_L2$scenario_rec <- ifelse(mlm_data_us_L2$scenario == 4, 1, 
                                      ifelse(mlm_data_us_L2$scenario == 1, 4,
                                             ifelse(mlm_data_us_L2$scenario == 2, 2, 3)))

#Recoding condition instead of using dummy variables to make visualization easier (coding human as reference cat).

mlm_data_us_L2$condition_rec <- ifelse(mlm_data_us_L2$condition == 2, 0, 
                                       ifelse(mlm_data_us_L2$condition == 1, 1, 2))

model_e1_us <- lme(value ~ factor(condition_rec)*factor(scenario_rec) + scale(EDU),
                   random =~ 1|ID, data=mlm_data_us_L2, control = lmeControl(opt = "optim"))

# Getting the estimated marginal means
emm1 = emmeans::emmeans(model_e1_us, ~ factor(condition_rec) * factor(scenario_rec))
plot(emm1)
emmip(model_e1_us, condition_rec ~ scenario_rec)
emm1

#AI as the reference category -- RQ1 (Figure 2)

model_e1_us_AI <- lme(value ~ Human + AI_Human + moderate + civil_reminder + present_facts +
                        Human*moderate + Human*civil_reminder + Human*present_facts + 
                        AI_Human*moderate + AI_Human*civil_reminder + AI_Human*present_facts + scale(EDU),
                      random =~ 1|ID, data=mlm_data_us_L2)

summary(model_e1_us_AI)
anova(model_e1_us_AI)

## Reestimated models with present facts as the reference category for Figure 2

###Multilevel Models for Figure 2

#Human as the reference category 

model_e1_us_facts <- lme(value ~ AI + AI_Human + moderate + civil_reminder + news +
                           AI*moderate + AI*civil_reminder + AI*news + 
                           AI_Human*moderate + AI_Human*civil_reminder + AI_Human*news + scale(EDU),
                         random =~ 1|ID, data=mlm_data_us_L2)
summary(model_e1_us_facts)
anova(model_e1_us_facts)


model_e1_us_AIfacts <- lme(value ~ Human + AI_Human + moderate + civil_reminder + news +
                             Human*moderate + Human*civil_reminder + Human*news + 
                             AI_Human*moderate + AI_Human*civil_reminder + AI_Human*news + EDU,
                           random =~ 1|ID, data=mlm_data_us_L2)

summary(model_e1_us_AIfacts)
anova(model_e1_us_AIfacts)


#############################################################################################
#                                                                                           #
#          2.3 Multilevel models for experiment-1 for pre-registered outcomes               #
#                     (Appendix 10. Experiment 1)                                           #
#                                                                                           #
#############################################################################################


####################################### FAIR ################################################################################################
### 3	Convert the wide data format to long data as MLM uses wide data format. Since, we have four scenarios we
###       consider each scenario as level. hence the data has four levels for experiment 1/

mlm_data_us_fair <- melt(mlm_data_us, id.vars = c("ID", "condition", "EDU") , measure.vars = c("FAIR", "CIVIL_FAIR", "FACTS_FAIR", "DEC_FAIR"))

mlm_data_us_fair$scenario <- recode(mlm_data_us_fair$variable, '"FAIR" = 1; "CIVIL_FAIR" = 2; "FACTS_FAIR" = 3; "DEC_FAIR" = 4;')

#Create dummy variables for conditions

mlm_data_us_fair$AI <- ifelse(mlm_data_us_fair$condition == 1, 1, 0)
mlm_data_us_fair$AI_Human <- ifelse(mlm_data_us_fair$condition == 3, 1, 0)
mlm_data_us_fair$Human <- ifelse(mlm_data_us_fair$condition == 2, 1, 0)
mlm_data_us_fair$moderate <- ifelse(mlm_data_us_fair$scenario == 1, 1, 0)
mlm_data_us_fair$civil_reminder <- ifelse(mlm_data_us_fair$scenario == 2, 1, 0)
mlm_data_us_fair$present_facts <- ifelse(mlm_data_us_fair$scenario == 3, 1, 0)
mlm_data_us_fair$news <- ifelse(mlm_data_us_fair$scenario == 4, 1, 0)

#Human as the reference category (Appendix 10. Experiment 1. Table 1. MLM results for the US)

model_e1_us_fair <- lme(value ~ AI + AI_Human + moderate + civil_reminder + present_facts +
                          AI*moderate + AI*civil_reminder + AI*present_facts + 
                          AI_Human*moderate + AI_Human*civil_reminder + AI_Human*present_facts + scale(EDU),
                        random =~ 1|ID, data=mlm_data_us_fair)
summary(model_e1_us_fair)
anova(model_e1_us_fair)

#Recoding scenario instead of using dummy variables to make visualization easier (coding news as reference cat).
mlm_data_us_fair$scenario_rec <- ifelse(mlm_data_us_fair$scenario == 4, 1, 
                                        ifelse(mlm_data_us_fair$scenario == 1, 4,
                                               ifelse(mlm_data_us_fair$scenario == 2, 2, 3)))

#Recoding condition instead of using dummy variables to make visualization easier (coding human as reference cat).

mlm_data_us_fair$condition_rec <- ifelse(mlm_data_us_fair$condition == 2, 0, 
                                         ifelse(mlm_data_us_fair$condition == 1, 1, 2))

#Model with recoded dummy variables

model_e1_us_fair <- lme(value ~ factor(condition_rec)*factor(scenario_rec) + scale(EDU),
                        random =~ 1|ID, data=mlm_data_us_fair, control = lmeControl(opt = "optim"))

# Getting the estimated marginal means (Appendix 10. Experiment 1. Table 1. Marginal means for the US)
emm1_fair = emmeans::emmeans(model_e1_us_fair, ~ factor(condition_rec) * factor(scenario_rec))
plot(emm1_fair)
emmip(model_e1_us_fair, condition_rec ~ scenario_rec)
emm1_fair

####################################### BIAS ################################################################################################
### 3	Convert the wide data format to long data as MLM uses wide data format. Since, we have four scenarios we
###       consider each scenario as level. hence the data has four levels for experiment 1

mlm_data_us_bias <- melt(mlm_data_us, id.vars = c("ID", "condition", "EDU") , measure.vars = c("BIASED", "CIVIL_BIASED", "FACTS_BIASED", "DEC_BIASED"))

mlm_data_us_bias$scenario <- recode(mlm_data_us_bias$variable, '"BIASED" = 1; "CIVIL_BIASED" = 2; "FACTS_BIASED" = 3; "DEC_BIASED" = 4;')

#Create dummy variables for conditions

mlm_data_us_bias$AI <- ifelse(mlm_data_us_bias$condition == 1, 1, 0)
mlm_data_us_bias$AI_Human <- ifelse(mlm_data_us_bias$condition == 3, 1, 0)
mlm_data_us_bias$Human <- ifelse(mlm_data_us_bias$condition == 2, 1, 0)
mlm_data_us_bias$moderate <- ifelse(mlm_data_us_bias$scenario == 1, 1, 0)
mlm_data_us_bias$civil_reminder <- ifelse(mlm_data_us_bias$scenario == 2, 1, 0)
mlm_data_us_bias$present_facts <- ifelse(mlm_data_us_bias$scenario == 3, 1, 0)
mlm_data_us_bias$news <- ifelse(mlm_data_us_bias$scenario == 4, 1, 0)

#Human as the reference category (Appendix 10. Experiment 1. Table 1. MLM results for the US)

model_e1_us_bias <- lme(value ~ AI + AI_Human + moderate + civil_reminder + present_facts +
                          AI*moderate + AI*civil_reminder + AI*present_facts + 
                          AI_Human*moderate + AI_Human*civil_reminder + AI_Human*present_facts + scale(EDU),
                        random =~ 1|ID, data=mlm_data_us_bias)
summary(model_e1_us_bias)
anova(model_e1_us_bias)

#Recoding scenario instead of using dummy variables to make visualization easier (coding news as reference cat).
mlm_data_us_bias$scenario_rec <- ifelse(mlm_data_us_bias$scenario == 4, 1, 
                                        ifelse(mlm_data_us_bias$scenario == 1, 4,
                                               ifelse(mlm_data_us_bias$scenario == 2, 2, 3)))

#Recoding condition instead of using dummy variables to make visualization easier (coding human as reference cat).

mlm_data_us_bias$condition_rec <- ifelse(mlm_data_us_bias$condition == 2, 0, 
                                         ifelse(mlm_data_us_bias$condition == 1, 1, 2))

#Model with recoded dummy variables

model_e1_us_bias <- lme(value ~ factor(condition_rec)*factor(scenario_rec) + scale(EDU),
                        random =~ 1|ID, data=mlm_data_us_bias, control = lmeControl(opt = "optim"))

# Getting the estimated marginal means (Appendix 10. Experiment 1. Table 1. Marginal means for the US)
emm1_bias = emmeans::emmeans(model_e1_us_bias, ~ factor(condition_rec) * factor(scenario_rec))
plot(emm1_bias)
emmip(model_e1_us_bias, condition_rec ~ scenario_rec)
emm1_bias

####################################### TRUST ################################################################################################
### 3	Convert the wide data format to long data as MLM uses wide data format. Since, we have four scenarios we
###       consider each scenario as level. hence the data has four levels for experiment 1

mlm_data_us_trust <- melt(mlm_data_us, id.vars = c("ID", "condition", "EDU") , measure.vars = c("TRUST", "CIVIL_TRUST", "FACTS_TRUST", "DEC_TRUST"))

mlm_data_us_trust$scenario <- recode(mlm_data_us_trust$variable, '"TRUST" = 1; "CIVIL_TRUST" = 2; "FACTS_TRUST" = 3; "DEC_TRUST" = 4;')

#Create dummy variables for conditions

mlm_data_us_trust$AI <- ifelse(mlm_data_us_trust$condition == 1, 1, 0)
mlm_data_us_trust$AI_Human <- ifelse(mlm_data_us_trust$condition == 3, 1, 0)
mlm_data_us_trust$Human <- ifelse(mlm_data_us_trust$condition == 2, 1, 0)
mlm_data_us_trust$moderate <- ifelse(mlm_data_us_trust$scenario == 1, 1, 0)
mlm_data_us_trust$civil_reminder <- ifelse(mlm_data_us_trust$scenario == 2, 1, 0)
mlm_data_us_trust$present_facts <- ifelse(mlm_data_us_trust$scenario == 3, 1, 0)
mlm_data_us_trust$news <- ifelse(mlm_data_us_trust$scenario == 4, 1, 0)

#Human as the reference category (Appendix 10. Experiment 1. Table 1. MLM results for the US)

model_e1_us_trust <- lme(value ~ AI + AI_Human + moderate + civil_reminder + present_facts +
                           AI*moderate + AI*civil_reminder + AI*present_facts + 
                           AI_Human*moderate + AI_Human*civil_reminder + AI_Human*present_facts + scale(EDU),
                         random =~ 1|ID, data=mlm_data_us_trust)
summary(model_e1_us_trust)
anova(model_e1_us_trust)

#Recoding scenario instead of using dummy variables to make visualization easier (coding news as reference cat).
mlm_data_us_trust$scenario_rec <- ifelse(mlm_data_us_trust$scenario == 4, 1, 
                                         ifelse(mlm_data_us_trust$scenario == 1, 4,
                                                ifelse(mlm_data_us_trust$scenario == 2, 2, 3)))

#Recoding condition instead of using dummy variables to make visualization easier (coding human as reference cat).

mlm_data_us_trust$condition_rec <- ifelse(mlm_data_us_trust$condition == 2, 0, 
                                          ifelse(mlm_data_us_trust$condition == 1, 1, 2))

#Model with recoded dummy variables

model_e1_us_trust <- lme(value ~ factor(condition_rec)*factor(scenario_rec) + scale(EDU),
                         random =~ 1|ID, data=mlm_data_us_trust, control = lmeControl(opt = "optim"))

# Getting the estimated marginal means (Appendix 10. Experiment 1. Table 1. Marginal means for the US)
emm1_trust = emmeans::emmeans(model_e1_us_trust, ~ factor(condition_rec) * factor(scenario_rec))
plot(emm1_trust)
emmip(model_e1_us_trust, condition_rec ~ scenario_rec)
emm1_trust

####################################### AGREE ################################################################################################
mlm_data_us_agree <- melt(mlm_data_us, id.vars = c("ID", "condition", "EDU") , measure.vars = c("AGREE", "CIVIL_AGREE", "FACTS_AGREE", "DEC_AGREE"))

mlm_data_us_agree$scenario <- recode(mlm_data_us_agree$variable, '"AGREE" = 1; "CIVIL_AGREE" = 2; "FACTS_AGREE" = 3; "DEC_AGREE" = 4;')

#Create dummy variables for conditions

mlm_data_us_agree$AI <- ifelse(mlm_data_us_agree$condition == 1, 1, 0)
mlm_data_us_agree$AI_Human <- ifelse(mlm_data_us_agree$condition == 3, 1, 0)
mlm_data_us_agree$Human <- ifelse(mlm_data_us_agree$condition == 2, 1, 0)
mlm_data_us_agree$moderate <- ifelse(mlm_data_us_agree$scenario == 1, 1, 0)
mlm_data_us_agree$civil_reminder <- ifelse(mlm_data_us_agree$scenario == 2, 1, 0)
mlm_data_us_agree$present_facts <- ifelse(mlm_data_us_agree$scenario == 3, 1, 0)
mlm_data_us_agree$news <- ifelse(mlm_data_us_agree$scenario == 4, 1, 0)

#Human as the reference category (Appendix 10. Experiment 1. Table 1. MLM results for the US)

model_e1_us_agree <- lme(value ~ AI + AI_Human + moderate + civil_reminder + present_facts +
                           AI*moderate + AI*civil_reminder + AI*present_facts + 
                           AI_Human*moderate + AI_Human*civil_reminder + AI_Human*present_facts + scale(EDU),
                         random =~ 1|ID, data=mlm_data_us_agree)
summary(model_e1_us_agree)
anova(model_e1_us_agree)

#Recoding scenario instead of using dummy variables to make visualization easier (coding news as reference cat).
mlm_data_us_agree$scenario_rec <- ifelse(mlm_data_us_agree$scenario == 4, 1, 
                                         ifelse(mlm_data_us_agree$scenario == 1, 4,
                                                ifelse(mlm_data_us_agree$scenario == 2, 2, 3)))

#Recoding condition instead of using dummy variables to make visualization easier (coding human as reference cat).

mlm_data_us_agree$condition_rec <- ifelse(mlm_data_us_agree$condition == 2, 0, 
                                          ifelse(mlm_data_us_agree$condition == 1, 1, 2))

#Model with recoded dummy variables

model_e1_us_agree <- lme(value ~ factor(condition_rec)*factor(scenario_rec) + scale(EDU),
                         random =~ 1|ID, data=mlm_data_us_agree, control = lmeControl(opt = "optim"))

# Getting the estimated marginal means (Appendix 10. Experiment 1. Table 1. Marginal means for the US)
emm1_agree = emmeans::emmeans(model_e1_us_agree, ~ factor(condition_rec) * factor(scenario_rec))
plot(emm1_agree)
emmip(model_e1_us_agree, condition_rec ~ scenario_rec)
emm1_agree

#############################################################################################
#                                                                                           #
#                   3  Multilevel models for experiment-2                                   #
#                                                                                           #
#                                                                                           #
#############################################################################################
########################################################################
#                                                                      #
#                   3.1 Creating variables of interest                 #
#                                                                      #
########################################################################

#### creating dummy variables (IVs) for AI, Humnan, and Mix conditions 
final_data$AI.dummy <- ifelse(final_data$condition %in% c(2,3),0,1)
final_data$Human.dummy <- ifelse(final_data$condition %in% c(1,3),0,1)
final_data$AI_Human.dummy <- ifelse(final_data$condition %in% c(1,2),0,1)
# Here AI condition is 1, Human condition is 2, and Mix condition is 3

final_data$counter_pro <- ifelse(final_data$pro %in% c(1),0,1)  



########################################################################
#                                                                      #
#                   3.2 Preparing data for MLM models                  #
#                                                                      #
########################################################################
### creating data backup
mlm_data_us <- final_data

### Creating VIF Function to check the multi-colinearity
vif.lme <- function (fit) {
  ## adapted from rms::vif
  v <- vcov(fit)
  nam <- names(fixef(fit))
  ## exclude intercepts
  ns <- sum(1 * (nam == "Intercept" | nam == "(Intercept)"))
  if (ns > 0) {
    v <- v[-(1:ns), -(1:ns), drop = FALSE]
    nam <- nam[-(1:ns)] }
  d <- diag(v)^0.5
  v <- diag(solve(v/(d %o% d)))
  names(v) <- nam
  v }

### Computing moderators
## the only moderator to be used in figures is mlm_data_us$avgppid; the rest in appendix tables

mlm_data_us$avgppid <- scale((mlm_data_us$PPID_1 + mlm_data_us$PPID_2 + mlm_data_us$PPID_3 +mlm_data_us$PPID_4)/4)
mlm_data_us$partystr <- scale(abs(mlm_data_us$PARTY7-4))
mlm_data_us$ideostr <- scale(abs(mlm_data_us$IDEO_1-5))
mlm_data_us$benefits_str = scale(abs(mlm_data_us$BEN_3_1 + mlm_data_us$BEN_3_2 + mlm_data_us$BEN_3_3)/3 - 6.33)
mlm_data_us$immig_str = scale(abs(mlm_data_us$IMM_3_1 + mlm_data_us$IMM_3_2 + mlm_data_us$IMM_3_3)/3 - 6.33)
mlm_data_us$econ_str = scale(abs(mlm_data_us$TRUMP_ECON_3_1 + mlm_data_us$TRUMP_ECON_3_2 +  mlm_data_us$TRUMP_ECON_3_3)/3 - 6.33)
mlm_data_us$counter_pro <- ifelse(mlm_data_us$pro ==1,0,1)

### 2 Load the data and assign the unique ID to participants. Generate the aggregated DV for all scenarios.

#creating datset backup
#assinging ID to each participant
mlm_data_us$ID <- 1:nrow(mlm_data_us)
mlm_data_us$agg_s1 <- (mlm_data_us$IMM_CRED + mlm_data_us$IMM_QUALITY + mlm_data_us$IMM_AGREE + mlm_data_us$IMM_BIAS)/4
mlm_data_us$agg_s2 <- (mlm_data_us$BEN_CRED + mlm_data_us$BEN_QUALITY + mlm_data_us$BEN_AGREE + mlm_data_us$BEN_BIAS)/4
mlm_data_us$agg_s3 <- (mlm_data_us$TE_CRED + mlm_data_us$TE_QUALITY + mlm_data_us$TE_AGREE + mlm_data_us$TE_BIAS)/4


### 3	Convert the wide data format to long data as MLM uses wide data format. Since, we have three issues we
###       consider each issue as level. hence the data has three levels for study 2.

mlm_data_us_L <- melt(mlm_data_us, id.vars = c("ID", "EDU", "counter_pro",
                                               "avgppid", "partystr", "ideostr", "benefits_str", "immig_str", "econ_str"), 
                      measure.vars = c("agg_s1", "agg_s2", "agg_s3"))
mlm_data_us_L2 <- melt(mlm_data_us, id.vars = c("ID", "condition", "EDU", "counter_pro",
                                                "avgppid", "partystr", "ideostr", "benefits_str", "immig_str", "econ_str"), 
                       measure.vars = c("agg_s1", "agg_s2", "agg_s3"))

library(car)
mlm_data_us_L2$scenario <- recode(mlm_data_us_L2$variable, '"agg_s1" = 1; "agg_s2" = 2;
"agg_s3" = 3;')

#Create dummy variables for conditions

mlm_data_us_L2$AI <- ifelse(mlm_data_us_L2$condition == 1, 1, 0)
mlm_data_us_L2$Human <- ifelse(mlm_data_us_L2$condition == 2, 1, 0)
mlm_data_us_L2$AI_Human <- ifelse(mlm_data_us_L2$condition == 3, 1, 0)
mlm_data_us_L2$immigration <- ifelse(mlm_data_us_L2$scenario == 1, 1, 0)
mlm_data_us_L2$ben <- ifelse(mlm_data_us_L2$scenario == 2, 1, 0)
mlm_data_us_L2$te <- ifelse(mlm_data_us_L2$scenario == 3, 1, 0)

########################################################################
#                                                                      #
#            3.3 MLM models main analysis                              #
#                                                                      #
########################################################################

#Multilevel baseline Model H3 -- is counter seen more negative than pro (Appendix 9 table 1)

model_e2_us_H3 <- lme(value ~ counter_pro + immigration + ben + scale(EDU) +
                        immigration*counter_pro + ben*counter_pro,
                      random =~ 1|ID, data=mlm_data_us_L2, control = lmeControl(opt = "optim"))
summary(model_e2_us_H3)
anova(model_e2_us_H3)
vif.lme(model_e2_us_H3)


#Multilevel baseline Model H4 -- ONLY ON THE COUNTER-ATT CONDITION TO AVOID three-way interactions and high values of 
# multicollinearity VIF = 8.15, although still below 10)

## AI as the reference (Appendix 9 table 2)
model_e2_us_counter <- lme(value ~ Human + AI_Human + immigration + ben + scale(EDU) +
                             Human*immigration + Human*ben +
                             AI_Human*immigration + AI_Human*ben,
                           random =~1|ID, data=mlm_data_us_L2[mlm_data_us_L2$counter_pro == 1,], control = lmeControl(opt = "optim"))
summary(model_e2_us_counter)
anova(model_e2_us_counter)
vif.lme(model_e2_us_counter)
car::vifmodel_e2_us_counter

# Mix as reference category
model_e2_us_mix <- lme(value ~ AI + Human + immigration + ben + scale(EDU) +
                         AI*immigration + AI*ben +
                         Human*immigration + Human*ben,
                       random =~1|ID, data=mlm_data_us_L2[mlm_data_us_L2$counter_pro == 1,], control = lmeControl(opt = "optim"))
summary(model_e2_us_mix)
anova(model_e2_us_mix)
vif.lme(model_e2_us_mix)
car::vif(model_e2_us_mix)

## EXPLORATORY ANALYSIS ## to test if people more open to pro/counter from ai/human/mixed

model_e2_us_test <- lme(value ~ Human + AI_Human + counter_pro + scale(EDU) +
                          Human*counter_pro + AI_Human*counter_pro,
                        random =~ 1|ID, data=mlm_data_us_L2, control = lmeControl(opt = "optim"))
summary(model_e2_us_test)
anova(model_e2_us_test)
vif.lme(model_e2_us_test)

model_e2_us_test_AI <- lme(value ~ AI + AI_Human + counter_pro + scale(EDU) +
                             AI*counter_pro + AI_Human*counter_pro,
                           random =~ 1|ID, data=mlm_data_us_L2, control = lmeControl(opt = "optim"))
summary(model_e2_us_test_AI)
anova(model_e2_us_test_AI)
vif.lme(model_e2_us_test_AI)


###  H 7 heterogeneous effects (Appendix 12. Table 1)
### for pol identity strength for figures
model_e2_us_avgppid <- lme(value ~ Human + AI_Human + immigration + ben + scale(EDU) +
                             Human*immigration + Human*ben +
                             AI_Human*immigration + AI_Human*ben +
                             avgppid + avgppid*Human + avgppid*AI_Human + 
                             avgppid*Human*ben + avgppid*AI_Human*ben + avgppid*Human*immigration + avgppid*AI_Human*immigration,
                           random =~ 1|ID, data=mlm_data_us_L2[mlm_data_us_L2$counter_pro == 1,], 
                           control = lmeControl(opt = "optim"))
summary(model_e2_us_avgppid)
anova(model_e2_us_avgppid)
vif.lme(model_e2_us_avgppid)

#### for pol identity strength human as reference 
model_e2_us_avgppid_humanref <- lme(value ~ AI + AI_Human + immigration + ben + scale(EDU) +
                                      AI*immigration + AI*ben +
                                      AI_Human*immigration + AI_Human*ben +
                                      avgppid + avgppid*AI + avgppid*AI_Human + 
                                      avgppid*AI*ben + avgppid*AI_Human*ben + avgppid*AI*immigration + avgppid*AI_Human*immigration,
                                    random =~ 1|ID, data=mlm_data_us_L2[mlm_data_us_L2$counter_pro == 1,], 
                                    control = lmeControl(opt = "optim"))
summary(model_e2_us_avgppid_humanref)
anova(model_e2_us_avgppid_humanref)
vif.lme(model_e2_us_avgppid_humanref)

## for tables in the other moderators
# party strength as moderator
model_e2_us_partystr <- lme(value ~ Human + AI_Human + immigration + ben + scale(EDU) +
                              partystr + partystr*Human + partystr*AI_Human,
                            random =~ 1|ID, data=mlm_data_us_L2[mlm_data_us_L2$counter_pro == 1,], 
                            control = lmeControl(opt = "optim"))
summary(model_e2_us_partystr)
anova(model_e2_us_partystr)
vif.lme(model_e2_us_partystr)

## Ideology strength as a moderator (Appendix 13 table 1)
model_e2_us_ideostr <- lme(value ~ Human + AI_Human + immigration + ben + scale(EDU) +
                             Human*immigration + Human*ben +
                             AI_Human*immigration + AI_Human*ben +
                             ideostr + ideostr*Human + ideostr*AI_Human +
                             ideostr*Human*ben + ideostr*AI_Human*ben + ideostr*Human*immigration + ideostr*AI_Human*immigration,
                           random =~ 1|ID, data=mlm_data_us_L2[mlm_data_us_L2$counter_pro == 1,], 
                           control = lmeControl(opt = "optim"))
summary(model_e2_us_ideostr)
anova(model_e2_us_ideostr)
vif.lme(model_e2_us_ideostr)

# H7 Heterogeneous effects issue attitudes (Interacting scenarios with issue att_strt, only for counter condition)
#(Appendix 13 table 2)
model_e2_us_att <- lme(value ~ Human + AI_Human + immigration + ben + scale(EDU) +
                         immig_str + benefits_str +econ_str +
                         Human*immigration + Human*ben +  AI_Human*immigration + AI_Human*ben +
                         Human*immig_str + AI_Human*immig_str +
                         Human*benefits_str + AI_Human*benefits_str + Human*econ_str + AI_Human*econ_str,
                       random =~ 1|ID, data=mlm_data_us_L2[mlm_data_us_L2$counter_pro == 1,], 
                       control = lmeControl(opt = "optim"))

summary(model_e2_us_att)
anova(model_e2_us_att)
vif.lme(model_e2_us_att)


########################################################################
#                                                                      #
#   3.4 MLM models per DV for the experiment2 (Appendix 11 table 1)    #
#                                                                      #
########################################################################

# creating table per DV for experiment 2 for US
mlm_data_us_e2_agree <- melt(mlm_data_us, id.vars = c("ID", "condition", "EDU", "counter_pro") , measure.vars = c("IMM_AGREE", "BEN_AGREE", "TE_AGREE"))
mlm_data_us_e2_agree$scenario <- recode(mlm_data_us_e2_agree$variable, '"IMM_AGREE" = 1; "BEN_AGREE" = 2; "TE_AGREE" = 3;')

mlm_data_us_e2_agree$AI <- ifelse(mlm_data_us_e2_agree$condition == 1, 1, 0)
mlm_data_us_e2_agree$Human <- ifelse(mlm_data_us_e2_agree$condition == 2, 1, 0)
mlm_data_us_e2_agree$AI_Human <- ifelse(mlm_data_us_e2_agree$condition == 3, 1, 0)
mlm_data_us_e2_agree$immigration <- ifelse(mlm_data_us_e2_agree$scenario == 1, 1, 0)
mlm_data_us_e2_agree$ben <- ifelse(mlm_data_us_e2_agree$scenario == 2, 1, 0)
mlm_data_us_e2_agree$te <- ifelse(mlm_data_us_e2_agree$scenario == 3, 1, 0)

model_e2_us_H4_agree <- lme(value ~ Human + AI_Human + immigration + ben + scale(EDU) +
                              Human*immigration + Human*ben +
                              AI_Human*immigration + AI_Human*ben,
                            random =~1|ID, data=mlm_data_us_e2_agree[mlm_data_us_e2_agree$counter_pro == 1,], control = lmeControl(opt = "optim"))
summary(model_e2_us_H4_agree)

mlm_data_us_e2_quality <- melt(mlm_data_us, id.vars = c("ID", "condition", "EDU", "counter_pro") , measure.vars = c("IMM_QUALITY", "BEN_QUALITY", "TE_QUALITY"))
mlm_data_us_e2_quality$scenario <- recode(mlm_data_us_e2_quality$variable, '"IMM_QUALITY" = 1; "BEN_QUALITY" = 2; "TE_QUALITY" = 3;')

mlm_data_us_e2_quality$AI <- ifelse(mlm_data_us_e2_quality$condition == 1, 1, 0)
mlm_data_us_e2_quality$Human <- ifelse(mlm_data_us_e2_quality$condition == 2, 1, 0)
mlm_data_us_e2_quality$AI_Human <- ifelse(mlm_data_us_e2_quality$condition == 3, 1, 0)
mlm_data_us_e2_quality$immigration <- ifelse(mlm_data_us_e2_quality$scenario == 1, 1, 0)
mlm_data_us_e2_quality$ben <- ifelse(mlm_data_us_e2_quality$scenario == 2, 1, 0)
mlm_data_us_e2_quality$te <- ifelse(mlm_data_us_e2_quality$scenario == 3, 1, 0)

model_e2_us_H4_quality <- lme(value ~ Human + AI_Human + immigration + ben + scale(EDU) +
                                Human*immigration + Human*ben +
                                AI_Human*immigration + AI_Human*ben,
                              random =~1|ID, data=mlm_data_us_e2_quality[mlm_data_us_e2_quality$counter_pro == 1,], control = lmeControl(opt = "optim"))
summary(model_e2_us_H4_quality)


mlm_data_us_e2_cred <- melt(mlm_data_us, id.vars = c("ID", "condition", "EDU", "counter_pro") , measure.vars = c("IMM_CRED", "BEN_CRED", "TE_CRED"))
mlm_data_us_e2_cred$scenario <- recode(mlm_data_us_e2_cred$variable, '"IMM_CRED" = 1; "BEN_CRED" = 2; "TE_CRED" = 3;')

mlm_data_us_e2_cred$AI <- ifelse(mlm_data_us_e2_cred$condition == 1, 1, 0)
mlm_data_us_e2_cred$Human <- ifelse(mlm_data_us_e2_cred$condition == 2, 1, 0)
mlm_data_us_e2_cred$AI_Human <- ifelse(mlm_data_us_e2_cred$condition == 3, 1, 0)
mlm_data_us_e2_cred$immigration <- ifelse(mlm_data_us_e2_cred$scenario == 1, 1, 0)
mlm_data_us_e2_cred$ben <- ifelse(mlm_data_us_e2_cred$scenario == 2, 1, 0)
mlm_data_us_e2_cred$te <- ifelse(mlm_data_us_e2_cred$scenario == 3, 1, 0)

model_e2_us_H4_cred <- lme(value ~ Human + AI_Human + immigration + ben + scale(EDU) +
                             Human*immigration + Human*ben +
                             AI_Human*immigration + AI_Human*ben,
                           random =~1|ID, data=mlm_data_us_e2_cred[mlm_data_us_e2_cred$counter_pro == 1,], control = lmeControl(opt = "optim"))
summary(model_e2_us_H4_cred)




mlm_data_us_e2_bias <- melt(mlm_data_us, id.vars = c("ID", "condition", "EDU", "counter_pro") , measure.vars = c("IMM_BIASED", "BEN_BIASED", "TE_BIASED"))
mlm_data_us_e2_bias$scenario <- recode(mlm_data_us_e2_bias$variable, '"IMM_BIASED" = 1; "BEN_BIASED" = 2; "TE_BIASED" = 3;')

mlm_data_us_e2_bias$AI <- ifelse(mlm_data_us_e2_bias$condition == 1, 1, 0)
mlm_data_us_e2_bias$Human <- ifelse(mlm_data_us_e2_bias$condition == 2, 1, 0)
mlm_data_us_e2_bias$AI_Human <- ifelse(mlm_data_us_e2_bias$condition == 3, 1, 0)
mlm_data_us_e2_bias$immigration <- ifelse(mlm_data_us_e2_bias$scenario == 1, 1, 0)
mlm_data_us_e2_bias$ben <- ifelse(mlm_data_us_e2_bias$scenario == 2, 1, 0)
mlm_data_us_e2_bias$te <- ifelse(mlm_data_us_e2_bias$scenario == 3, 1, 0)
model_e2_us_H4_bias <- lme(value ~ Human + AI_Human + immigration + ben + scale(EDU) +
                             Human*immigration + Human*ben +
                             AI_Human*immigration + AI_Human*ben,
                           random =~1|ID, data=mlm_data_us_e2_bias[mlm_data_us_e2_bias$counter_pro == 1,], control = lmeControl(opt = "optim"))
summary(model_e2_us_H4_bias)


########################################################################
#                                                                      #
#   3.5 Composite Marginal means for the experiment2                   #
#                                                                      #
########################################################################


mlm_data_us_L2$scenario <- recode(mlm_data_us_L2$variable, '"agg_s1" = 1; "agg_s2" = 2;
"agg_s3" = 3;')

#Recoding scenario instead of using dummy variables to make visualization easier (coding trump economy as reference cat).
mlm_data_us_L2$scenario_rec <- ifelse(mlm_data_us_L2$scenario == 3, 1, ifelse(mlm_data_us_L2$scenario == 1, 3, 2))


#Model for H3, H4 and H5 using factor variables instead of dummy (results are the same, only variable names are less clear in output)
#(Appendix-9 table 3)
model_e2_us_emm <- lme(value ~ factor(condition) * factor(scenario_rec) * factor(counter_pro) + scale(EDU),
                       random =~ 1|ID, data=mlm_data_us_L2, control = lmeControl(opt = "optim"))

### Condition 1 = AI, Condition 2 = Human, Condition 3 = Mixed
### Issue 1 = Trump economy, Issue 2 = Benefits, Issue 3 = Immigration
### Counter_pro 0 = pro-attitudinal, Counter_pro 1 = counter-attitudinal

summary(model_e2_us_emm)
anova(model_e2_us_emm)
vif.lme(model_e2_us_emm)

### Getting the estimated marginal means
emm1 = emmeans::emmeans(model_e2_us_emm, ~ factor(condition) * factor(scenario_rec) * factor(counter_pro))
emm1


#### H7 Heterogeneous effects

### Party ID (Appendix 12 table 2)

model_e2_us_avgppid <- lme(value ~ factor(condition) * factor(scenario_rec) * avgppid + scale(EDU),
                           random =~ 1|ID, data=mlm_data_us_L2[mlm_data_us_L2$counter_pro == 1,], 
                           control = lmeControl(opt = "optim"))
summary(model_e2_us_avgppid)
anova(model_e2_us_avgppid)
vif.lme(model_e2_us_avgppid)

### Getting the estimated marginal means.
#NOTE: Because we have a continupus moderator, we need the argument cov.reduce = range
emm2 = emmeans::emmeans(model_e2_us_avgppid, ~ factor(condition) * factor(scenario_rec) * avgppid, cov.reduce = range)
emm2

### Party Strength

model_e2_us_partystr <- lme(value ~ factor(condition) * factor(scenario_rec) * partystr + scale(EDU),
                            random =~ 1|ID, data=mlm_data_us_L2[mlm_data_us_L2$counter_pro == 1,], 
                            control = lmeControl(opt = "optim"))
summary(model_e2_us_partystr)
anova(model_e2_us_partystr)
vif.lme(model_e2_us_partystr)

# Estimated Marginal Means
emm3 = emmeans::emmeans(model_e2_us_partystr, ~ factor(condition) * factor(scenario_rec) * partystr, cov.reduce = range)
emm3
